So here we’re going to do some simpler model selection things with the data using GLMs instead of GAMs
Sets up the functions that will be used later
rad2deg <- function(rad) {(rad * 180) / (pi)}
deg2rad <- function(deg) {(deg * pi) / (180)}
round_any <- function(x, accuracy, f=round){f(x/ accuracy) * accuracy}
ang_mean <- function(x){rad2deg(atan(mean(sin(deg2rad(x)))/mean(cos(deg2rad(x)))))}
fold_angle_0_360_to_0_180 <- function(x){abs(abs(x-180)-180)}
Reads in the data and alters it as needed
comp_data <- read.csv("Fish_Comp_Values.csv")
comp_data <- na.omit(comp_data)
comp_data <- comp_data %>% mutate(Flow = ifelse(Flow == "0", "Flow 0", "Flow 2")) %>%
mutate(Ablation = ifelse(Ablation == "N", "No Ablation", "Ablated")) %>%
mutate(Darkness = ifelse(Darkness == "N", "Light", "Dark")) %>%
mutate(Heading_Diff = abs(Heading_Diff)) %>%
filter(Distance <= 4) %>%
filter(Speed_Diff <= 6) %>%
mutate(quarter_heading_diff = abs(abs(Heading_Diff-90)-90)) %>%
mutate(Is_Aligned = ifelse(Heading_Diff < 30, 1, 0)) %>%
mutate(Is_Reversed = ifelse(Heading_Diff > 150, 1, 0)) %>%
mutate(Angle = fold_angle_0_360_to_0_180(Angle))
sum_comp_data <- comp_data %>% mutate(X_Distance = round_any(X_Distance,0.25),
Y_Distance = round_any(abs(Y_Distance),0.25)) %>%
group_by(Flow,Ablation,Darkness,X_Distance,Y_Distance) %>%
summarise(Speed_Diff = mean(Speed_Diff),
Heading_Diff = ang_mean(Heading_Diff),
Sync = mean(Sync),
Quarter_Heading_Diff = mean(quarter_heading_diff),
Is_Aligned = mean(Is_Aligned),
Is_Reversed = mean(Is_Reversed))
`summarise()` has grouped output by 'Flow', 'Ablation', 'Darkness', 'X_Distance'. You can override using the `.groups` argument.
#Basic Stats
speed_anova <- glm(Speed_Diff ~ Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness, data = comp_data)
Anova(speed_anova)
Analysis of Deviance Table (Type II tests)
Response: Speed_Diff
LR Chisq Df Pr(>Chisq)
Flow 20.1114 1 7.306e-06 ***
Ablation 24.4652 1 7.566e-07 ***
Darkness 13.7065 1 0.0002137 ***
Flow:Ablation 0.0032 1 0.9545865
Flow:Darkness 1.5440 1 0.2140271
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
heading_anova <- glm(quarter_heading_diff ~ Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness, data = comp_data)
Anova(heading_anova)
Analysis of Deviance Table (Type II tests)
Response: quarter_heading_diff
LR Chisq Df Pr(>Chisq)
Flow 0.39037 1 0.5321
Ablation 0.11881 1 0.7303
Darkness 0.44816 1 0.5032
Flow:Ablation 0.44387 1 0.5053
Flow:Darkness 0.32170 1 0.5706
sync_anova <- glm(Sync ~ Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness, data = comp_data)
Anova(sync_anova)
Analysis of Deviance Table (Type II tests)
Response: Sync
LR Chisq Df Pr(>Chisq)
Flow 2.0144 1 0.1558116
Ablation 14.0452 1 0.0001785 ***
Darkness 0.7801 1 0.3771189
Flow:Ablation 9.3099 1 0.0022792 **
Flow:Darkness 7.8318 1 0.0051335 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dist_anova <- glm(Distance ~ Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness, data = comp_data)
Anova(dist_anova)
Analysis of Deviance Table (Type II tests)
Response: Distance
LR Chisq Df Pr(>Chisq)
Flow 0.358 1 0.5496
Ablation 0.068 1 0.7944
Darkness 107.389 1 <2e-16 ***
Flow:Ablation 1.970 1 0.1605
Flow:Darkness 0.032 1 0.8577
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
ggplot(comp_data, aes(x = Flow, y = Speed_Diff, fill = paste(Ablation,Darkness,sep="-")))+
geom_boxplot(outlier.shape = NA) +
ylim(0,3) +
theme_light()
ggplot(comp_data, aes(x = Flow, y = quarter_heading_diff, fill = paste(Ablation,Darkness,sep="-")))+
geom_boxplot(outlier.shape = NA) +
ylim(0,90) +
theme_light()
ggplot(comp_data, aes(x = Flow, y = Sync, fill = paste(Ablation,Darkness,sep="-")))+
geom_boxplot(outlier.shape = NA) +
ylim(0,1) +
theme_light()
ggplot(comp_data, aes(x = Flow, y = Distance, fill = paste(Ablation,Darkness,sep="-")))+
geom_boxplot(outlier.shape = NA) +
ylim(0,4) +
theme_light()
ggplot(comp_data, aes(x = Flow, y = Angle, fill = paste(Ablation,Darkness,sep="-")))+
geom_boxplot(outlier.shape = NA) +
ylim(0,180) +
theme_light()
Now we add in distance and angle
speed_anova_da <- glm(Speed_Diff ~ poly(Distance,4)*Angle*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)
Anova(speed_anova_da)
Analysis of Deviance Table (Type II tests)
Response: Speed_Diff
LR Chisq Df Pr(>Chisq)
poly(Distance, 4) 115.177 4 < 2.2e-16 ***
Angle 5.946 1 0.0147510 *
Flow 27.012 1 2.022e-07 ***
Ablation 27.929 1 1.258e-07 ***
Darkness 6.729 1 0.0094872 **
poly(Distance, 4):Angle 1.964 4 0.7423022
Flow:Ablation 0.006 1 0.9399196
Flow:Darkness 0.300 1 0.5840184
poly(Distance, 4):Flow 1.872 4 0.7592434
poly(Distance, 4):Ablation 18.003 4 0.0012324 **
poly(Distance, 4):Darkness 20.232 4 0.0004494 ***
Angle:Flow 5.300 1 0.0213266 *
Angle:Ablation 1.855 1 0.1732289
Angle:Darkness 0.750 1 0.3863939
poly(Distance, 4):Flow:Ablation 24.820 4 5.468e-05 ***
poly(Distance, 4):Flow:Darkness 8.276 4 0.0819669 .
Angle:Flow:Ablation 0.000 1 0.9967220
Angle:Flow:Darkness 0.068 1 0.7937129
poly(Distance, 4):Angle:Flow 14.614 4 0.0055717 **
poly(Distance, 4):Angle:Ablation 12.121 4 0.0164745 *
poly(Distance, 4):Angle:Darkness 3.761 4 0.4393154
poly(Distance, 4):Angle:Flow:Ablation 26.352 4 2.687e-05 ***
poly(Distance, 4):Angle:Flow:Darkness 9.819 4 0.0435844 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
heading_anova_da <- glm(quarter_heading_diff ~ Distance*Angle*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)
Anova(heading_anova_da)
Analysis of Deviance Table (Type II tests)
Response: quarter_heading_diff
LR Chisq Df Pr(>Chisq)
Distance 2.5860 1 0.10781
Angle 0.0237 1 0.87761
Flow 0.5212 1 0.47035
Ablation 0.1588 1 0.69030
Darkness 0.3356 1 0.56238
Distance:Angle 2.1602 1 0.14162
Flow:Ablation 0.5623 1 0.45332
Flow:Darkness 0.1771 1 0.67388
Distance:Flow 0.0729 1 0.78715
Distance:Ablation 1.2093 1 0.27148
Distance:Darkness 0.1028 1 0.74850
Angle:Flow 0.0598 1 0.80688
Angle:Ablation 0.9390 1 0.33254
Angle:Darkness 0.9958 1 0.31832
Distance:Flow:Ablation 0.2560 1 0.61286
Distance:Flow:Darkness 1.4412 1 0.22995
Angle:Flow:Ablation 0.2370 1 0.62640
Angle:Flow:Darkness 0.0116 1 0.91409
Distance:Angle:Flow 0.1515 1 0.69713
Distance:Angle:Ablation 4.7465 1 0.02936 *
Distance:Angle:Darkness 0.7361 1 0.39092
Distance:Angle:Flow:Ablation 1.3680 1 0.24216
Distance:Angle:Flow:Darkness 6.3862 1 0.01150 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
sync_anova_da <- glm(Sync ~ Distance*Angle*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)
Anova(sync_anova_da)
Analysis of Deviance Table (Type II tests)
Response: Sync
LR Chisq Df Pr(>Chisq)
Distance 31.817 1 1.694e-08 ***
Angle 0.666 1 0.4143894
Flow 1.399 1 0.2368323
Ablation 13.757 1 0.0002081 ***
Darkness 0.029 1 0.8647383
Distance:Angle 3.098 1 0.0783799 .
Flow:Ablation 8.824 1 0.0029737 **
Flow:Darkness 7.255 1 0.0070689 **
Distance:Flow 1.241 1 0.2652318
Distance:Ablation 0.036 1 0.8487728
Distance:Darkness 1.186 1 0.2761968
Angle:Flow 2.525 1 0.1120840
Angle:Ablation 5.473 1 0.0193101 *
Angle:Darkness 1.156 1 0.2822980
Distance:Flow:Ablation 7.958 1 0.0047885 **
Distance:Flow:Darkness 2.776 1 0.0956635 .
Angle:Flow:Ablation 1.892 1 0.1689323
Angle:Flow:Darkness 0.068 1 0.7941314
Distance:Angle:Flow 0.201 1 0.6536996
Distance:Angle:Ablation 11.861 1 0.0005733 ***
Distance:Angle:Darkness 0.150 1 0.6983841
Distance:Angle:Flow:Ablation 0.871 1 0.3506822
Distance:Angle:Flow:Darkness 7.252 1 0.0070829 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
ggplot(comp_data, aes(x = Distance, y = Speed_Diff, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
geom_smooth(method = "glm", formula = y ~ poly(x,4))+
facet_wrap(~ Flow) +
theme_light()
ggplot(comp_data, aes(x = Distance, y = quarter_heading_diff, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
geom_smooth(method = "glm", formula = y ~ poly(x,4))+
facet_wrap(~ Flow) +
theme_light()
ggplot(comp_data, aes(x = Distance, y = Sync, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
geom_smooth(method = "glm", formula = y ~ poly(x,4))+
facet_wrap(~ Flow) +
theme_light()
ggplot(comp_data, aes(x = Angle, y = Speed_Diff, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
geom_smooth(method = "glm", formula = y ~ poly(x,4))+
facet_wrap(~ Flow) +
theme_light()
ggplot(comp_data, aes(x = Angle, y = quarter_heading_diff, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
geom_smooth(method = "glm", formula = y ~ poly(x,4))+
facet_wrap(~ Flow) +
theme_light()
ggplot(comp_data, aes(x = Angle, y = Sync, color = paste(Ablation,Darkness,sep="-"), fill = paste(Ablation,Darkness,sep="-")))+
geom_smooth(method = "glm", formula = y ~ poly(x,4))+
facet_wrap(~ Flow) +
theme_light()
This creates a dataframe that is used for making model predictions and graphing them
d <- seq(from = 0, to = 4, by = 0.1)
a <- seq(from = 0, to = 180, by = 10)
flows <- c("Flow 0", "Flow 2")
ablation <- c("No Ablation", "Ablated")
dark <- c("Light","Dark")
predict_df_da <- expand.grid(Distance = d, Angle = a, Flow = flows, Ablation = ablation, Darkness = dark)
predict_df_da <- predict_df_da %>% mutate(X_Distance = Distance*(cos(deg2rad(Angle))), Y_Distance = Distance*(sin(deg2rad(Angle))))
predict_df_da <- predict_df_da %>% filter(!(Ablation == "Ablated" & Darkness == 'Dark'))
predict_df_da <- na.omit(predict_df_da)
So I want to make a model that tries to include both distance and angle. What I really need is to figure out what power works best for angle and distance, and then include that. The issue is stepwise selection can just remove one of them entirely. So I’m going to use cross validation to select the power for both of them, and then use that. I will also calculate the AICs to help in this decision making process
set.seed(7)
comp_data_speed_model <- comp_data %>% select(c(Flow,Darkness,Ablation,Angle,Distance,Speed_Diff))
max_poly <- 10
speed_dist_cv_error_10 <- rep(0,max_poly)
speed_dist_AIC_10 <- rep(0,max_poly)
for (i in 1:max_poly){
speed_dist_fit <- glm(Speed_Diff ~ poly(Distance,i)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
data = comp_data_speed_model)
speed_dist_cv_error_10[i] <- cv.glm(comp_data_speed_model, speed_dist_fit, K = 10)$delta[1]
speed_dist_AIC_10[i] <- AIC(speed_dist_fit)
}
speed_angle_cv_error_10 <- rep(0,max_poly)
speed_angle_AIC_10 <- rep(0,max_poly)
for (i in 1:max_poly){
speed_angle_fit <- glm(Speed_Diff ~ poly(Angle,i)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
data = comp_data_speed_model)
speed_angle_cv_error_10[i] <- cv.glm(comp_data_speed_model, speed_angle_fit, K = 10)$delta[1]
speed_angle_AIC_10[i] <- AIC(speed_angle_fit)
}
speed_poly_plot_df <- data.frame(c(seq(max_poly),seq(max_poly)),
c(speed_dist_cv_error_10,speed_angle_cv_error_10),
c(speed_dist_AIC_10,speed_angle_AIC_10),
c(rep("Distance",max_poly),rep("Angle",max_poly)))
colnames(speed_poly_plot_df) <- c("Degree","Error","AIC","Predictor")
speed_poly_plot_df <- speed_poly_plot_df %>% group_by(Predictor) %>%
mutate(minError = min(Error),minAIC = min(AIC)) %>%
ungroup() %>%
mutate(isMinEror = ifelse(Error == minError,3,1),isMinAIC = ifelse(AIC == minAIC,3,1))
ggplot(speed_poly_plot_df, aes(x = Degree, y = Error, color = Predictor))+
geom_point(size = speed_poly_plot_df$isMinEror)+
geom_line()+
theme_light()+
facet_wrap(~ Predictor, scales = "free") +
scale_size(guide = "none")
ggplot(speed_poly_plot_df, aes(x = Degree, y = AIC, color = Predictor))+
geom_point(size = speed_poly_plot_df$isMinAIC)+
geom_line()+
theme_light()+
facet_wrap(~ Predictor, scales = "free") +
scale_size(guide = "none")
Now let’s use those polynomials to create the actual model
While they not be the lowest, the big improvements happen at 4 for Distance, and 2 for Angle
speed_diff_glm <- glm(Speed_Diff ~ poly(Distance,4)*poly(Angle,2)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
data = comp_data_speed_model)
summary(speed_diff_glm)
Call:
glm(formula = Speed_Diff ~ poly(Distance, 4) * poly(Angle, 2) *
(Ablation + Flow + Darkness + Ablation:Flow + Darkness:Flow),
data = comp_data_speed_model)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.2116 -0.3067 -0.0859 0.1910 3.4546
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.374e-01 3.163e-02 26.476 < 2e-16 ***
poly(Distance, 4)1 5.475e+00 3.118e+00 1.756 0.079188 .
poly(Distance, 4)2 6.853e+00 3.187e+00 2.150 0.031583 *
poly(Distance, 4)3 1.147e+00 2.850e+00 0.402 0.687379
poly(Distance, 4)4 7.173e+00 2.413e+00 2.973 0.002962 **
poly(Angle, 2)1 4.186e+00 2.370e+00 1.766 0.077374 .
poly(Angle, 2)2 -2.217e+00 2.752e+00 -0.806 0.420391
AblationNo Ablation 7.305e-02 2.163e-02 3.377 0.000736 ***
FlowFlow 2 5.948e-02 4.749e-02 1.253 0.210403
DarknessLight -4.160e-02 2.755e-02 -1.510 0.131132
poly(Distance, 4)1:poly(Angle, 2)1 4.232e+02 2.309e+02 1.832 0.066938 .
poly(Distance, 4)2:poly(Angle, 2)1 3.148e+02 2.369e+02 1.329 0.183886
poly(Distance, 4)3:poly(Angle, 2)1 3.883e+02 2.266e+02 1.713 0.086702 .
poly(Distance, 4)4:poly(Angle, 2)1 -9.528e+01 2.035e+02 -0.468 0.639653
poly(Distance, 4)1:poly(Angle, 2)2 -1.237e+02 3.028e+02 -0.409 0.682821
poly(Distance, 4)2:poly(Angle, 2)2 2.186e+02 3.088e+02 0.708 0.479134
poly(Distance, 4)3:poly(Angle, 2)2 3.175e+02 2.579e+02 1.231 0.218397
poly(Distance, 4)4:poly(Angle, 2)2 8.375e+01 2.127e+02 0.394 0.693829
AblationNo Ablation:FlowFlow 2 1.443e-02 3.753e-02 0.384 0.700689
FlowFlow 2:DarknessLight 4.637e-03 4.152e-02 0.112 0.911073
poly(Distance, 4)1:AblationNo Ablation -5.285e+00 2.275e+00 -2.323 0.020220 *
poly(Distance, 4)2:AblationNo Ablation -2.449e+00 2.303e+00 -1.063 0.287643
poly(Distance, 4)3:AblationNo Ablation -1.821e+00 2.002e+00 -0.910 0.362897
poly(Distance, 4)4:AblationNo Ablation -3.468e+00 1.604e+00 -2.162 0.030663 *
poly(Distance, 4)1:FlowFlow 2 -7.224e+00 4.631e+00 -1.560 0.118819
poly(Distance, 4)2:FlowFlow 2 -1.121e+01 4.809e+00 -2.330 0.019838 *
poly(Distance, 4)3:FlowFlow 2 -4.884e-01 4.382e+00 -0.111 0.911259
poly(Distance, 4)4:FlowFlow 2 -4.360e-01 3.740e+00 -0.117 0.907201
poly(Distance, 4)1:DarknessLight 2.302e+00 2.659e+00 0.866 0.386686
poly(Distance, 4)2:DarknessLight -5.323e+00 2.769e+00 -1.922 0.054593 .
poly(Distance, 4)3:DarknessLight -1.956e+00 2.514e+00 -0.778 0.436673
poly(Distance, 4)4:DarknessLight -4.435e+00 2.144e+00 -2.069 0.038619 *
poly(Angle, 2)1:AblationNo Ablation -3.353e+00 1.624e+00 -2.065 0.039006 *
poly(Angle, 2)2:AblationNo Ablation 4.392e+00 1.976e+00 2.223 0.026278 *
poly(Angle, 2)1:FlowFlow 2 -1.248e+00 3.654e+00 -0.342 0.732653
poly(Angle, 2)2:FlowFlow 2 8.118e+00 4.250e+00 1.910 0.056166 .
poly(Angle, 2)1:DarknessLight -1.067e+00 2.076e+00 -0.514 0.607152
poly(Angle, 2)2:DarknessLight -6.070e-01 2.352e+00 -0.258 0.796389
poly(Distance, 4)1:AblationNo Ablation:FlowFlow 2 1.004e+01 3.834e+00 2.617 0.008880 **
poly(Distance, 4)2:AblationNo Ablation:FlowFlow 2 4.709e+00 3.926e+00 1.199 0.230398
poly(Distance, 4)3:AblationNo Ablation:FlowFlow 2 1.692e-02 3.463e+00 0.005 0.996103
poly(Distance, 4)4:AblationNo Ablation:FlowFlow 2 -2.326e-01 2.877e+00 -0.081 0.935552
poly(Distance, 4)1:FlowFlow 2:DarknessLight 4.597e+00 3.940e+00 1.167 0.243451
poly(Distance, 4)2:FlowFlow 2:DarknessLight 1.025e+01 4.104e+00 2.498 0.012525 *
poly(Distance, 4)3:FlowFlow 2:DarknessLight -5.024e-01 3.813e+00 -0.132 0.895185
poly(Distance, 4)4:FlowFlow 2:DarknessLight 5.413e-01 3.249e+00 0.167 0.867682
poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2 1.924e+00 2.904e+00 0.663 0.507641
poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2 -6.275e+00 3.508e+00 -1.789 0.073740 .
poly(Angle, 2)1:FlowFlow 2:DarknessLight 8.569e-01 3.212e+00 0.267 0.789669
poly(Angle, 2)2:FlowFlow 2:DarknessLight 1.495e+00 3.659e+00 0.409 0.682838
poly(Distance, 4)1:poly(Angle, 2)1:AblationNo Ablation -3.168e+02 1.646e+02 -1.925 0.054294 .
poly(Distance, 4)2:poly(Angle, 2)1:AblationNo Ablation -4.343e+02 1.636e+02 -2.655 0.007962 **
poly(Distance, 4)3:poly(Angle, 2)1:AblationNo Ablation -9.630e+01 1.515e+02 -0.636 0.525092
poly(Distance, 4)4:poly(Angle, 2)1:AblationNo Ablation 5.509e+01 1.300e+02 0.424 0.671750
poly(Distance, 4)1:poly(Angle, 2)2:AblationNo Ablation 1.351e+02 2.312e+02 0.584 0.559017
poly(Distance, 4)2:poly(Angle, 2)2:AblationNo Ablation -5.181e+01 2.345e+02 -0.221 0.825118
poly(Distance, 4)3:poly(Angle, 2)2:AblationNo Ablation -1.465e+02 1.912e+02 -0.766 0.443742
poly(Distance, 4)4:poly(Angle, 2)2:AblationNo Ablation 1.054e+02 1.488e+02 0.708 0.478770
poly(Distance, 4)1:poly(Angle, 2)1:FlowFlow 2 -7.079e+02 3.478e+02 -2.035 0.041884 *
poly(Distance, 4)2:poly(Angle, 2)1:FlowFlow 2 -9.543e+01 3.550e+02 -0.269 0.788075
poly(Distance, 4)3:poly(Angle, 2)1:FlowFlow 2 -5.019e+02 3.380e+02 -1.485 0.137620
poly(Distance, 4)4:poly(Angle, 2)1:FlowFlow 2 1.022e+03 3.057e+02 3.344 0.000832 ***
poly(Distance, 4)1:poly(Angle, 2)2:FlowFlow 2 3.442e+02 4.573e+02 0.753 0.451688
poly(Distance, 4)2:poly(Angle, 2)2:FlowFlow 2 -4.818e+02 4.694e+02 -1.026 0.304771
poly(Distance, 4)3:poly(Angle, 2)2:FlowFlow 2 -3.571e+02 4.016e+02 -0.889 0.373974
poly(Distance, 4)4:poly(Angle, 2)2:FlowFlow 2 -3.991e+02 3.306e+02 -1.207 0.227334
poly(Distance, 4)1:poly(Angle, 2)1:DarknessLight -7.976e+01 1.999e+02 -0.399 0.689922
poly(Distance, 4)2:poly(Angle, 2)1:DarknessLight 1.168e+02 2.100e+02 0.556 0.578190
poly(Distance, 4)3:poly(Angle, 2)1:DarknessLight -2.198e+02 2.025e+02 -1.086 0.277713
poly(Distance, 4)4:poly(Angle, 2)1:DarknessLight 1.097e+02 1.827e+02 0.601 0.548090
poly(Distance, 4)1:poly(Angle, 2)2:DarknessLight 1.008e+02 2.534e+02 0.398 0.690692
poly(Distance, 4)2:poly(Angle, 2)2:DarknessLight -3.290e+02 2.654e+02 -1.239 0.215268
poly(Distance, 4)3:poly(Angle, 2)2:DarknessLight -1.985e+02 2.247e+02 -0.883 0.377124
poly(Distance, 4)4:poly(Angle, 2)2:DarknessLight -2.237e+02 1.881e+02 -1.189 0.234394
poly(Distance, 4)1:poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2 2.673e+02 2.841e+02 0.941 0.346800
poly(Distance, 4)2:poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2 6.864e+01 2.836e+02 0.242 0.808736
poly(Distance, 4)3:poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2 -6.060e+01 2.623e+02 -0.231 0.817321
poly(Distance, 4)4:poly(Angle, 2)1:AblationNo Ablation:FlowFlow 2 -9.584e+02 2.302e+02 -4.163 3.18e-05 ***
poly(Distance, 4)1:poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2 -3.361e+02 3.930e+02 -0.855 0.392475
poly(Distance, 4)2:poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2 1.646e+02 3.993e+02 0.412 0.680267
poly(Distance, 4)3:poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2 7.080e+01 3.353e+02 0.211 0.832742
poly(Distance, 4)4:poly(Angle, 2)2:AblationNo Ablation:FlowFlow 2 1.890e+02 2.615e+02 0.723 0.469682
poly(Distance, 4)1:poly(Angle, 2)1:FlowFlow 2:DarknessLight 3.462e+02 2.970e+02 1.166 0.243717
poly(Distance, 4)2:poly(Angle, 2)1:FlowFlow 2:DarknessLight -1.392e+02 3.029e+02 -0.460 0.645818
poly(Distance, 4)3:poly(Angle, 2)1:FlowFlow 2:DarknessLight 4.288e+02 2.946e+02 1.456 0.145562
poly(Distance, 4)4:poly(Angle, 2)1:FlowFlow 2:DarknessLight -6.988e+02 2.675e+02 -2.612 0.009022 **
poly(Distance, 4)1:poly(Angle, 2)2:FlowFlow 2:DarknessLight 1.313e+02 3.807e+02 0.345 0.730104
poly(Distance, 4)2:poly(Angle, 2)2:FlowFlow 2:DarknessLight 4.487e+02 3.929e+02 1.142 0.253424
poly(Distance, 4)3:poly(Angle, 2)2:FlowFlow 2:DarknessLight 2.905e+02 3.426e+02 0.848 0.396498
poly(Distance, 4)4:poly(Angle, 2)2:FlowFlow 2:DarknessLight 5.870e+02 2.864e+02 2.050 0.040424 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.2194725)
Null deviance: 1416.0 on 6050 degrees of freedom
Residual deviance: 1308.3 on 5961 degrees of freedom
AIC: 8086.8
Number of Fisher Scoring iterations: 2
Now let’s make some predictions
speed_pred <- predict_df_da %>% mutate(Speed_Diff = predict(speed_diff_glm,predict_df_da))
comp_data <- comp_data %>% mutate(Round_Dist = as.factor(round_any(Distance,1)), Round_Angle = as.factor(round_any(Angle,30)))
ggplot()+
geom_boxplot(data = comp_data, aes(x = Round_Dist, y = Speed_Diff))+
facet_wrap(~ Flow + Ablation + Darkness) +
theme_light()
round_dist_aov <- aov(Speed_Diff ~ Round_Dist*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)
Anova(round_dist_aov)
Anova Table (Type II tests)
Response: Speed_Diff
Sum Sq Df F value Pr(>F)
Round_Dist 21.92 4 24.3669 < 2.2e-16 ***
Flow 5.30 1 23.5595 1.242e-06 ***
Ablation 5.98 1 26.5765 2.614e-07 ***
Darkness 1.68 1 7.4726 0.006283 **
Flow:Ablation 0.00 1 0.0057 0.939866
Flow:Darkness 0.10 1 0.4649 0.495347
Round_Dist:Flow 0.96 4 1.0642 0.372465
Round_Dist:Ablation 2.19 4 2.4329 0.045321 *
Round_Dist:Darkness 2.99 4 3.3206 0.010029 *
Round_Dist:Flow:Ablation 5.31 4 5.9036 9.730e-05 ***
Round_Dist:Flow:Darkness 1.94 4 2.1605 0.070840 .
Residuals 1354.08 6021
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
ggplot()+
geom_boxplot(data = comp_data, aes(x = Round_Angle, y = Speed_Diff))+
facet_wrap(~ Flow + Ablation + Darkness) +
theme_light()
round_angle_aov <- aov(Speed_Diff ~ Round_Angle*(Flow + Ablation + Darkness + Flow:Ablation + Flow:Darkness), data = comp_data)
Anova(round_angle_aov)
Anova Table (Type II tests)
Response: Speed_Diff
Sum Sq Df F value Pr(>F)
Round_Angle 3.25 6 2.3765 0.0269962 *
Flow 4.90 1 21.5029 3.608e-06 ***
Ablation 6.21 1 27.2924 1.808e-07 ***
Darkness 2.99 1 13.1415 0.0002912 ***
Flow:Ablation 0.05 1 0.2310 0.6307756
Flow:Darkness 0.36 1 1.5723 0.2099197
Round_Angle:Flow 7.43 6 5.4377 1.280e-05 ***
Round_Angle:Ablation 2.02 6 1.4765 0.1818843
Round_Angle:Darkness 1.43 6 1.0461 0.3930796
Round_Angle:Flow:Ablation 5.99 6 4.3881 0.0001975 ***
Round_Angle:Flow:Darkness 0.98 6 0.7154 0.6372082
Residuals 1368.14 6009
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ggplot()+
# geom_point(data = comp_data %>% filter(Speed_Diff <= 2), aes(x = Distance, y = Speed_Diff), alpha = 0.1)+
# geom_smooth(data = speed_pred, aes(x = Distance, y = Speed_Diff))+
# facet_wrap(~ Flow + Ablation + Darkness) +
# theme_light()
# ggplot()+
# geom_density_2d_filled(data = comp_data %>% filter(Speed_Diff <= 2), aes(x = Distance, y = Speed_Diff), contour_var = "ndensity")+
# geom_smooth(data = speed_pred, aes(x = Distance, y = Speed_Diff, color = "red"))+
# facet_wrap(~ Flow + Ablation + Darkness) +
# theme_light()
# ggplot()+
# geom_point(data = comp_data %>% filter(Speed_Diff <= 2), aes(x = Angle, y = Speed_Diff), alpha = 0.1)+
# geom_smooth(data = speed_pred, aes(x = Angle, y = Speed_Diff))+
# facet_wrap(~ Flow + Ablation + Darkness) +
# theme_light()
# ggplot()+
# geom_density_2d_filled(data = comp_data %>% filter(Speed_Diff <= 2), aes(x = Angle, y = Speed_Diff), contour_var = "ndensity")+
# geom_smooth(data = speed_pred, aes(x = Angle, y = Speed_Diff, color = "red"))+
# facet_wrap(~ Flow + Ablation + Darkness) +
# theme_light()
ggplot(data = comp_data, aes(x = Angle, y = Distance, z = Speed_Diff))+
stat_summary_2d() +
facet_wrap(~ Flow + Ablation + Darkness) +
scale_fill_viridis(direction = -1) +
theme_light()
You know let’s just try stepwise as well. There’s jsut way too many thing in that possible model
speed_m_all <- glm(Speed_Diff ~ Distance*Angle*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow) +
I(Distance^2)*I(Angle^2)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow)+
I(Distance^3)*I(Angle^3)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow)+
I(Distance^4)*I(Angle^4)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
data = comp_data_speed_model)
speed_m_none <- glm(Speed_Diff ~ 1, data = comp_data_speed_model)
speed_m_both <- step(speed_m_none, direction = "both", scope = formula(speed_m_all), trace = F)
summary(speed_m_both)
Call:
glm(formula = Speed_Diff ~ Flow + Distance + I(Angle^4) + Ablation +
Angle + I(Angle^3) + I(Angle^2) + Flow:I(Angle^4) + Distance:Ablation +
Flow:Ablation + Ablation:Angle + Flow:Angle + Flow:Distance +
Flow:Distance:Ablation, data = comp_data_speed_model)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.8986 -0.3680 -0.1121 0.2432 3.7254
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.535e-01 6.916e-02 2.220 0.026451 *
FlowFlow 2 2.876e-01 6.243e-02 4.608 4.16e-06 ***
Distance 1.190e-01 1.528e-02 7.791 7.77e-15 ***
I(Angle^4) -4.649e-09 1.644e-09 -2.828 0.004706 **
AblationNo Ablation 3.053e-01 5.171e-02 5.905 3.73e-09 ***
Angle 8.898e-03 3.795e-03 2.345 0.019083 *
I(Angle^3) 1.682e-06 6.132e-07 2.743 0.006113 **
I(Angle^2) -1.965e-04 7.740e-05 -2.539 0.011145 *
FlowFlow 2:I(Angle^4) 4.946e-10 1.137e-10 4.350 1.38e-05 ***
Distance:AblationNo Ablation -9.326e-02 1.966e-02 -4.743 2.15e-06 ***
FlowFlow 2:AblationNo Ablation -1.349e-01 6.483e-02 -2.081 0.037497 *
AblationNo Ablation:Angle -1.082e-03 3.179e-04 -3.404 0.000668 ***
FlowFlow 2:Angle -2.065e-03 6.409e-04 -3.222 0.001280 **
FlowFlow 2:Distance -4.181e-02 2.399e-02 -1.743 0.081447 .
FlowFlow 2:Distance:AblationNo Ablation 1.065e-01 3.036e-02 3.508 0.000456 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.2679219)
Null deviance: 1716.9 on 6103 degrees of freedom
Residual deviance: 1631.4 on 6089 degrees of freedom
AIC: 9300
Number of Fisher Scoring iterations: 2
Now let’s try the same thing for sync values (Rayleigh’s R)
set.seed(7)
comp_data_sync_model <- comp_data %>% select(c(Flow,Darkness,Ablation,Angle,Distance,Sync))
max_poly <- 10
sync_dist_cv_error_10 <- rep(0,max_poly)
sybc_dist_AIC_10 <- rep(0,max_poly)
for (i in 1:max_poly){
sync_dist_fit <- glm(Sync ~ poly(Distance,i)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
data = comp_data_sync_model)
sync_dist_cv_error_10[i] <- cv.glm(comp_data_sync_model, sync_dist_fit, K = 10)$delta[1]
sybc_dist_AIC_10[i] <- AIC(sync_dist_fit)
}
sync_angle_cv_error_10 <- rep(0,max_poly)
sync_angle_AIC_10 <- rep(0,max_poly)
for (i in 1:max_poly){
sync_angle_fit <- glm(Sync ~ poly(Angle,i)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
data = comp_data_sync_model)
sync_angle_cv_error_10[i] <- cv.glm(comp_data_sync_model, sync_angle_fit, K = 10)$delta[1]
sync_angle_AIC_10[i] <- AIC(sync_angle_fit)
}
sync_poly_plot_df <- data.frame(c(seq(max_poly),seq(max_poly)),
c(speed_dist_cv_error_10,speed_angle_cv_error_10),
c(speed_dist_AIC_10,speed_angle_AIC_10),
c(rep("Distance",max_poly),rep("Angle",max_poly)))
colnames(sync_poly_plot_df) <- c("Degree","Error","AIC","Predictor")
sync_poly_plot_df <- sync_poly_plot_df %>% group_by(Predictor) %>%
mutate(minError = min(Error),minAIC = min(AIC)) %>%
ungroup() %>%
mutate(isMinEror = ifelse(Error == minError,3,1),isMinAIC = ifelse(AIC == minAIC,3,1))
ggplot(sync_poly_plot_df, aes(x = Degree, y = Error, color = Predictor))+
geom_point(size = sync_poly_plot_df$isMinEror)+
geom_line()+
theme_light()+
facet_wrap(~ Predictor, scales = "free") +
scale_size(guide = "none")
ggplot(sync_poly_plot_df, aes(x = Degree, y = AIC, color = Predictor))+
geom_point(size = sync_poly_plot_df$isMinAIC)+
geom_line()+
theme_light()+
facet_wrap(~ Predictor, scales = "free") +
scale_size(guide = "none")
Now let’s use those polynomials to create the actual model
While they not be the lowest, the big improvements happen at 2 for Distance, and 2 for Angle
sync_diff_glm <- glm(Sync ~ I(Distance^2)*I(Angle^2)*(Ablation+Flow+Darkness+Ablation:Flow+Darkness:Flow),
data = comp_data_sync_model)
summary(sync_diff_glm)
Call:
glm(formula = Sync ~ I(Distance^2) * I(Angle^2) * (Ablation +
Flow + Darkness + Ablation:Flow + Darkness:Flow), data = comp_data_sync_model)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.55625 -0.19623 -0.01846 0.19226 0.51838
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.374e-01 3.651e-02 11.983 < 2e-16 ***
I(Distance^2) 8.645e-03 6.445e-03 1.341 0.179875
I(Angle^2) 5.721e-06 2.576e-06 2.221 0.026404 *
AblationNo Ablation 8.300e-02 2.232e-02 3.718 0.000202 ***
FlowFlow 2 2.099e-01 5.287e-02 3.971 7.24e-05 ***
DarknessLight 5.984e-02 3.305e-02 1.811 0.070222 .
I(Distance^2):I(Angle^2) -7.040e-07 5.514e-07 -1.277 0.201717
AblationNo Ablation:FlowFlow 2 -1.648e-01 3.754e-02 -4.389 1.16e-05 ***
FlowFlow 2:DarknessLight -9.831e-02 4.694e-02 -2.094 0.036263 *
I(Distance^2):AblationNo Ablation -1.454e-03 4.438e-03 -0.328 0.743149
I(Distance^2):FlowFlow 2 -3.488e-02 9.445e-03 -3.692 0.000224 ***
I(Distance^2):DarknessLight -1.244e-02 5.622e-03 -2.214 0.026893 *
I(Angle^2):AblationNo Ablation -3.000e-06 1.686e-06 -1.780 0.075200 .
I(Angle^2):FlowFlow 2 -1.480e-05 3.804e-06 -3.891 0.000101 ***
I(Angle^2):DarknessLight -3.767e-06 2.288e-06 -1.647 0.099688 .
I(Distance^2):AblationNo Ablation:FlowFlow 2 2.164e-02 7.423e-03 2.915 0.003567 **
I(Distance^2):FlowFlow 2:DarknessLight 2.284e-02 8.053e-03 2.836 0.004585 **
I(Angle^2):AblationNo Ablation:FlowFlow 2 8.795e-06 2.705e-06 3.251 0.001157 **
I(Angle^2):FlowFlow 2:DarknessLight 7.684e-06 3.377e-06 2.275 0.022920 *
I(Distance^2):I(Angle^2):AblationNo Ablation -4.831e-08 4.062e-07 -0.119 0.905336
I(Distance^2):I(Angle^2):FlowFlow 2 3.123e-06 8.161e-07 3.826 0.000131 ***
I(Distance^2):I(Angle^2):DarknessLight 7.798e-07 4.726e-07 1.650 0.098985 .
I(Distance^2):I(Angle^2):AblationNo Ablation:FlowFlow 2 -1.543e-06 6.555e-07 -2.354 0.018600 *
I(Distance^2):I(Angle^2):FlowFlow 2:DarknessLight -2.089e-06 7.046e-07 -2.965 0.003042 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.06028771)
Null deviance: 371.77 on 6103 degrees of freedom
Residual deviance: 366.55 on 6080 degrees of freedom
AIC: 204.5
Number of Fisher Scoring iterations: 2
Now let’s make some predictions
sync_pred <- predict_df_da %>% mutate(Sync = predict(sync_diff_glm,predict_df_da))
ggplot(comp_data, aes(x = Distance, y = Sync))+
geom_point(alpha = 0.1)+
geom_smooth(method = lm, formula = y ~ poly(x, 2)) +
facet_wrap(~ Flow + Ablation + Darkness) +
theme_light()
ggplot(comp_data, aes(x = Angle, y = Sync))+
geom_point(alpha = 0.1)+
geom_smooth(method = glm, formula = y ~ poly(x, 2))+
facet_wrap(~ Flow + Ablation + Darkness) +
theme_light()
ggplot(sync_pred, aes(x = Distance, y = Sync))+
geom_smooth()+
scale_fill_viridis() +
facet_wrap(~ Flow + Ablation + Darkness) +
theme_light()
ggplot(sync_pred, aes(x = Angle, y = Sync))+
geom_smooth()+
scale_fill_viridis() +
facet_wrap(~ Flow + Ablation + Darkness) +
theme_light()